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Abstract 

In modelling sediment compaction and mineral 
reactions, the rheological behaviour of sediments 
is typically considered as poroelastic or purely 
viscous. In fact, compaction due to pressure 
solution and mechanical processes in porous media 
is far more complicated. A generalised model 
of viscoelastic compaction and the smectite to 
illite mineral reaction in hydrocarbon basins is 
presented. A one-step dehydration model of the 
mineral reaction is assumed. The obtained non- 
linear governing equations are solved numerically 
and different combinations of physical parameters 
are used to simulate realistic situations in typical 
sedimentary basins. Comparison of numerical 
simulations with real data has shown very good 
agreement with respect to both the porosity profile 
and the mineral reaction. 

Key words: Compaction, mineral reaction, vis- 
coclasticity, sedimentary basins, pressure solution. 
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1 Introduction 

When a well-bore is being drilled for oil explo- 
ration, drilling mud is used in the hole to main- 
tain its integrity and safety. The mud density 
must be sufficient to prevent collapse of the hole, 
but not so high that hydrofracturing of the sur- 
rounding rock occurs. Both these effects depend 
on the pore fluid pressure in the rock, and drilling 
problems occur in regions where abnormal pore 
pressure occurs. Such overpressuring can substan- 
tially affect oil-drilling rates and even cause serious 
blowouts during drilling. Therefore, an industrially 
important objective is to predict overpressuring be- 
fore drilling and to identify its precursors during 
drilling. Another related objective is to predict 
reservoir quality and hydrocarbon migration. An 
essential step to achieve such objectives is the sci- 
entific understanding of overpressuring mechanisms 
and the evolutionary history of post-depositional 
sediments such as shales. 

The extent of compaction is strongly influenced 
by burial history and the lithology of sediments. 
The freshly deposited loosely packed sediments 
tend to evolve, like an open system, towards a 
closely packed grain framework during the initial 
stages of burial compaction and this is accom- 
plished by the processes of grain slippage, rota- 
tion, bending and brittle fracturing. Such reori- 
entation processes are collectively referred to as 
mechanical compaction, and generally take place 
in the first 1-2 km of burial. After this ini- 
tial porosity loss, further porosity reduction is ac- 
complished by processes of chemical compaction 
such as pressure solution (Rutter, 1983; Tada and 
Siever, 1989). Overpressuring may develop from 
a non- equilibrium compaction environment (Hed- 
berg, 1936; Gibson et al, 1967; Audet and Fowler, 
1992; Connolly, 1997; Connolly and Podladchikov, 
1998; Fowler and Yang, 1998). A rapidly accu- 
mulating basin is unable to expel pore fluids suf- 
ficiently rapidly due to the weight of overburden 
rock. The development of overpressuring retards 
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compaction, resulting in higher porosity, perme- 
ability and low thermal conductivity than are nor- 
mal for a given depth. These changes alter the 
structural and stratigraphic features of sedimen- 
tary units and provide the potential for hydrocar- 
bon migration. Rieke and Chilingarian (1974) gave 
a detailed review on this topic. Audet and Fowler 
(1992) also reviewed compaction briefly and Con- 
nolly and Podladchikov (2000) provided a more re- 
cent comprehensive review. 

Pressure solution has been considered as an im- 
portant process in deformation and porosity change 
during compaction in sedimentary rocks (Bathurst, 
1971;Tada and Sicvcr, 1989). Pressure solution 
refers to a process by which grains dissolve at in- 
tergranular contacts under non-hydrostatic stress 
and rcprecipitate in pore spaces, thus resulting in 
compaction. The solubility of minerals increases 
with increasing effective stress at grain contacts. 
Pressure dissolution at grain contacts is therefore 
a compactional response of the sediment during 
burial in an attempt to increase the grain contact 
area so as to distribute the effective stress over a 
larger surface. There has been renewed interest 
on pressure solution mechanisms because of its im- 
portant role in the diagenesis of sedimentary rocks 
and its relation to rock deformation. Two com- 
mon types of pressure solutions occur in nature. 
The intergranular pressure solution was first stud- 
ied by Sorby (1863) based while the stylolitiza- 
tion was first described by Fuchs (1894) followed 
by many authors such as Stockdale (1922), Dun- 
nington (1954) and Heald (1955). Tada and Siever 
(1989) gave a comprehensive review on pressure 
solution and some recent literature can be found 
in recent papers (Birchwood and Turcotte, 1994; 
Schneider, 1996; Connolly and Podladchikov, 1998; 
Fowler and Yang, 1999; Yang, 2000a). Pressure 
solution process is typically assumed to viscous 
(Weyl, 1959; Rutter, 1983; Birchwood and Tur- 
cotte, 1994; Fowler and Yang, 1999) and it is usu- 
ally referred to as viscous compaction, viscous creep 
or pressure solution creep. Its rheological constitu- 
tive relation (or compaction relation) is often writ- 
ten as a relationship between effective stress and 
strain rate. However, the more natural and real- 
istic rheology shall be viscoelastic (Connolly and 
Podladchikov, 1998; Revil, 1999; Yang, 2000c) and 
this is the approach we will use in the present pa- 
per. 

Mineral reactions are observed worldwide in sedi- 
mentary hydrocarbon basins. The close spatial cor- 
relations between smectite disappearance and illite 
formation imply the existence of the smectite-illite 



reaction. This transformation is one of the most 

important in clastic rocks. The reaction has re- 
ceived much attention but the nature of both the 
illite/smectite mixed-layer and the reaction mech- 
anism arc still under discussion (Eberl and Hower 
1976; Velde and Vasseur, 1992; Abercrombie et al., 
1994). Measured rates of this process in the lab- 
oratory (Eberl and Hower 1976) suggest that at 
elevated temperature, the transformation proceeds 
very fast from the geological point of view. On 
the other hand, observations suggest that this re- 
action is initiated relatively suddenly at a temper- 
ature of 90*^^7, but then takes place gradually over 
a depth of several himdrcd meters, which suggests 
a time scale of the order of a million years. The 
transformation can be considered a simple dehy- 
dration reaction. In fact, the mechanism of such a 
mineral reaction is rather more complicated and is 
not completely understood. The mineral reaction 
may consist of dissolution of montmorillonite in free 
pore water and subsequent precipitation of silica 
as illite. The following dissolution-precipitation re- 
actions have been proposed by Yang (2000b) and 
Fowler (2000) 
Smectite dissolution 

M^(smectite) ^ [Si'^-] + n[H20], (1) 

K-feldspar dissolution 

[K-feldspar] [K+^] + [AlO:^^-] + [SiO\], (2) 

Elite precipitation 

[K+^] + [AlO^^-] + [Si^-] r(illite) + [SiOl], 

(3) 

Quartz dissolution and precipitation 

[SiO!^] <H- [quartz], (4) 

where s, I denote solid and liquid phase. [—Si—] 
is an aqueous silica combination in such forms as 
[-{Si4)Oio{OH)2]. [A102^-] is only a general no- 
tation of the combination such as [Al{OH)^]. The 
reaction rates in the above dissolution and repre- 
cipitation mechanism may be quite different at dif- 
ferent steps because dissolution is usually enhanced 
by pressure solution due to the increased pressure 
at grain contacts, while precipipation at free pore 
space is less enhanced by pressure (Rimstidt and 
Barnes, 1980). In the limiting case, the above four 
step reactions can be approximately considered as 
a one-step (first order) dehydration process. The 
first-order dehydration model is a good approxima- 
tion in the sense of describing the extent of progress 
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of the overall smectite-to-illite transformation with- 
out much concern for its detailed geochemical fea- 
tures. Therefore, we represent it schematically as 

smectite ■ [-ff20]-^illite + nif2o[-ff20] (free water), 

(5) 

which was suggested in the earlier works (Eberl and 
Hower, 1976; Velde and Vasseur, 1992; Abercromie 
et al., 1994). This one-step model was studied 
briefly by Yang (2000b) and Fowler (2000). 

Despite the importance of compaction and min- 
eral reactions for geological problems, the literature 
dealing with quantitative modelling is not exten- 
sive. Although qualitative features have been re- 
ceived much attention in the literature, the mech- 
anismd are still poorly understood. A full under- 
standing of the mechanism requires an interdisci- 
plinary study involving soil mechanics, geochem- 
istry, geophysics and geology (Hedberg, 1936; Rieke 
and Chilingarian, 1974; Tada and Siever, 1989; 
Fowler and Yang, 1998; Connolly and Podlad- 
chikov,2000). This paper provides a unified com- 
paction relation in a form of a visco-poroelastic re- 
lation of Maxwell type. The nonlinear partial dif- 
ferential equations are then solved numerically and 
compared with real borehole log data. 



2 MATHEMATICAL 
MODEL 

For the convenience of investigating the effect of 
sediment compaction and mineral reactions, the 
sediment is considered cLS clS Qi compressible porous 
matrix. Combining mass conservation of the pore 
fluid with Darcy's law leads to model equations 
of the general type. Consider a matrix consist- 
ing of four interdispersed media: free pore wa- 
ter and three solid phase clay minerals, namely, 
quartz, water-rich montmorillonite and dehydrated 
illite. Let the volume fractions of the respec- 
tive media (montmorillonite, illite, quartz, water) 
be ^m, <j)i, <j)c, 0, so that 4)m + (t>i + 4>c + 4> = We 
assume that all the solids move with the same aver- 
aged velocity m", while; the pore water has velocity 
v} . The rate at which montmorillonite is trans- 
formed is denoted by r^, the rate of illite forma- 
tion is n , and the rate at which water released from 
the mineral reaction is = — r, due to mass 
conservation. 
Conservation of mass 



dt 



V ■ 



-^ + V-(</>u') = r^, 

+ + (l>i + (l)m = '^, 

Conservation of Energy 

d 



(7) 

(8) 
(9) 



dt 



+V • {[psCs{l - (^)u« + Pici(l>vi^]T} = V • {Kth^T), 

(10) 

Darcy's law 

<^(u'-U^) = --(Vp' + p;3j), (11) 

Force balance 
V-(T^-V[p']-P.9j = 0, p={l-^)ps + (t>Pu (12) 

where the momentum equation has been written in 
the similar form as Fowler (1990) and Audet and 
Fowler (1992). This present form is also equivalent 
to that given by (McKenzie, 1984) except for the 
notation difference, u' and are the velocities 
of fluid and solid matrix, k and p are the matrix 
permeability and the liquid viscosity, pi and ps are 
the densities of fluid and solid matrix, <r* is the 
effective stress, j is the unit vector pointing verti- 
cally upwards, p' is the pore pressure, and g is the 
gravitational acceleration. Kth is thermal conduc- 
tivity, and T is temperature, and q are specific 
heat of the solid and liquid, respectively. In for- 
mulating the mathematical model, the heat release 
from the mineral reaction has been neglected since 
it is usually very small. In addition, a rheological 
constitutive relation, which is often refered as com- 
paction law, is needed to complete this model and 
this is described in the following section 2.1 in more 
detail. 

Combing equation force balance and Darcy's law, 

we have 



</)(u'-u«) 



[V.a«-(p,-pO(l-«, (13) 



dt 



+ V • [</.„u^] 



(6) 



which is a reduced form of Darcy's law and j is a 
unit vector pointing upward. 

2.1 Poroelasticity, Pressure Solution 
and Viscoelastic Compaction 

Nonlinear compaction models have been formu- 
lated in two ways. The early and most common 
models assumed an elastic or poroelastic rheology. 
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and the compaction relation is of the Athy's type 
Pe = Pe{4>) (Gibson, England & Hussey, 1967; 
Smith, 1971; Sharp, 1976; Wangen, 1992; Audet 
and Fowler, 1992; Fowler and Yang, 1998). More 
recent models assumes a viscous rheology with a 
compaction relation of the form = — ^(</))V.u* 
where ^ is the bulk viscosity (Augevine and Tur- 
cottc. 1983; Birchwood and Turcottc, 1994; Fowler 
and Yang, 1999). The poroelastic models are valid 
for the mechanical compaction while the viscous 
models arc used to simulate the chemical com- 
paction such as pressure solution. We can thus 
generalise the above relations to a viscoelastic com- 
paction law of the 1-D Maxwell type 



1 Dpe 



1 



D__d_ 

dz ^ K{(j)) Dt eW*" Dt~dt^''d^.' 

(14) 

and 



K(0) 



(15) 



where Cg is a constant or compaction index de- 
scribing the degree of compaction of sediments, 
is the value of viscosity at </) = 0o (the initial 
porosity). Equivalently, we would anticipate a vis- 
coelastic rheology for the medium, involving mate- 
rial derivatives of tensors, and some care is needed 
to ensure that the resulting model is frame indif- 
ferent. That is to say, the rheological relation of 
stress-strain should be invariant under the coordi- 
nate transformation. This is not always guaranteed 
due to the complexity of the rheological relations 
(Bird, Armstrong & Hassager 1977). Fortunately, 
for one-dimensional flow, which is always irrota- 
tional, the equation is invariant and all the different 
equations in corotional and codeformational frames 
degenerate into the same form (Yang, 2000c). In 
the one-dimensional case we will discuss below, we 
can take this for granted. 



3 Non-dimensionalization 

For a 1-D basin, the physical domain is from the 
basin basement {z = 0) to the basin top (ocean 
floor) and the basin thickness varies with time t. 
The length-scale d is a typical length which will be 
defined 

P = i r-i> (16) 

[Ps- Pi)gd 

where G = 1 + \s a, constant describing sed- 
iment properties and Pe is the effective pressure. 



The choice of scaling d is in such a way that the di- 

mcnsionlcss pressure p = 0(1). Assigning the typ- 
ical thermal gradient to be 7, we also require that 
e = (T - To)/(7d) = 0(1). Meanwhile, we scale 
z with rf, with rhs, time t with d/ihs, perme- 
ability k with ko, and heat conductivity Kth with 
Kq where rhs,ko and Kq are set to typically ob- 
served values. By writing k{(f)) = kak*, z = dz*, 
Kth = KqK, and dropping the asterisks, we 
thus have 



dt 



-7^e^®</.„, , (17) 
7^(l - a)e'5«</.„, (18) 

= an5e^®(j)m, (19) 
= 1; (20) 
^{[a(l -</.)+ 0]e} + ^{[a(l - ct>)u^ + ct>u']Q) 



d(t>i ^ d{(t>iu'') 
dt dz 



dt dz 

4) + (f>i + (t>c 



dz dz 
dp 



(21) 



0(u'-«^) = AM<^)[-^-(l -</,)], (22) 



du' 
dz 



C,cj> Dp (</»/(/.o)" 

r-. Pi 



where 



(1 - .^)2 Dt 

ko{Ps - Pi)g 



(23) 



A = 



A = 



piciihsd' 
k2d 



5 = 



{Ps - pi)9(P 



Mm. 



Eg 



a = 



PsCs 



(24) 



(25) 



pi Pici 

and k^ and fco are the reaction rate and permeabil- 
ity at the basin top, respectively. Cg is a constant 
coefflcient or compaction index in the known func- 
tion K{(l)) in equation (14). is the molar wa- 
ter released from one molar smectite; Myj and 
are the molar weights of water and smectite, d is 
the typical length scale in the basin; Ea is the ac- 
tivation energy of the mineral reaction; Tq is the 
temperature at the ocean floor and R is the gas 
constant. In addition, ((/)o/^)" expresses the de- 
pendence of the bulk viscosity ^ on porosity, which 
has been determined empirically to be n = 1.3 for 
uncemented sand-like granular media and n = 2 for 
a wide range of cemented rocks (Paterson, 1995). 
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The nondimensional parameter A is essentially 

the ratio of hydroconductivity to the sedimenta- 
tion rate and thus characterises the relative time 
scales of fluid escape and sedimentation so that 
A <C 1 is the slow compaction limit while A ^ 1 is 
the fast compaction limit (Fowler and ang, 1998). 
Similarly, A is the ratio of the thermal conduc- 
tivity to scdimcntion rate, and thus controls the 
progress of the temperature evolution. The large 
A case (A ^ 1) so that a nearly equilibrium state 
is reached. The parameter TZ is the relative time 
scales of mineral reaction to sedimentation, and P 
is the reduced activation energy. However, the fac- 
tor TZcxp{pQ) always appears as a combination in 
the above model equations, and can be rewritten 
as 

7^exp(/3e) =exp[/3(e-ee)] and = ^\n^, 

(26) 

where the new parameter 9c, which replaces TZ, 

is a dimcnsionlcss critical temperature (with refer- 
ence to the surface temperature). In the following 
discussions, we will see that the smectite to illite 
reaction virtually takes place in a region called the 
reaction window, at a depth of ~ 0c, with its thick- 
ness controlled by /?. 

In the above derivation, we have used the re- 
quirement of degenerating to the poroelastic case 
p = \n.{(f)Q / (j)) — (00 — (t^) when S — >• oo. The consti- 
tutive relation for permeability k{(j)) is nonlinear, 
and its typical form (Smith, 1971) is 



where m{t) is the nondimensional sedimentation 

rate. As the thermal conductivity K{(f)) varies with 
(Nield,1991), a constant heat flux will leads to a 
varying temperature gradient. However, we will as- 
sume K = const (in fitting the borehole data) be- 
cause of its weak dependence on porosity and gives 
a good approximation in reality (Schneider, 1996). 

It is useful to estimate these parameters by 
using values taken from observations. By us- 
ing the typical values of pi ~ lO^kgm"^, ps ~ 



-3 



2.5 X lO'^kgm 

lO-^Nsm^ ^ 

300 m Ma" ^ 

c, - 500J Kg-i K-^ 
-1 j^-i 



-10-20 m2,/x 



1 X 10-11 m s-S g 



ci 



ko ~ 10-15 

1 X 10^1 N s m-2, TTis ~ 
10m s-2, 
4200JKg-i K-i, 
Ko - 1.5 Ww-i if-i. To ~ 280 K, G ~ 1, d ~ 
1000m, M„ ~ 450, Ea ~ SOkJmorSCs ~ 1.0; 
then A « 0.01 - 1000, A = 0.1-400, S w 1, /3 = 2.3 
and G = 3.0. An initial porosity of (po = 0.4 for 
pore water and initial volume fractions of 0.3 for 
montmorillonite, 0.0 for illite and 0.3 for quartz 
are used in the following computations. It is worth 
pointing out that the large variations of A and A 
are due to the fact that the permeability fco varies 
greatly over a wide range where the smaller values 
of ko correspond to a lower permeability and larger 
values correpond to higher peameability. This wide 
range can thus similate different types of sediments 
in different sedimentation environments (Rieke and 
ChiUngarian, 1974). 

4 Numerical Method 



(27) 



where m = 8 is a typical value for shales (Smith, 
1971; Audet and Fowler, 1992). However, the 
present model can use any value in the range from 
m = 3toTO = 11 in the numerical simulations. Ini- 
tially, the basin thickness is zero or h{t = 0) = 0. 
During the whole evolutionary process, new sedi- 
ment deposits at a constant rate rhs with a constant 
initial porosity (po- The basement rock is assumed 
to be impcamcablc or = = at z = while a 
constant heat flux {-KdQ/dz) is introduced at the 
basement rock. The top boundary is at constant 
pressure which can be taken to be zero. Thus, the 
boundary conditions become 



u' = 0, u' = 0, K— 

oz 



-1, at z = Q, (28) 



h = rn{t) + u^ at z = h{t). 



(29) 



In order to solve the highly coupled non-linear 
equations, an implicit numerical difference method 
is used. Substituting the expression for darcy flow 
into the other ciquations, the essential equations de- 
scribing for porosities (p and (pm become the stan- 
dard non- linear parabolic form (Meek and Norbury, 
1982) 

<Pt = F{z,t,cP)cP,,+g{z,t,cP,cP^). (30) 

The first step gives 0"+!/^ as a solution of the fol- 
lowing equation 



At 



(^)F(.„r+V2,^.)52^^+l/2 
+g{Zi,t^+^/^<P^,j-J,4>7), 



(31) 



where S^cpi = (0^+1 - 2(pi + (pi_i) and S^^cpi = 
(l/2)(0j+i — (pi-i). The second stage gives 0"'*'^ 
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as a solution of the following equation 

+g(.„r+i/2^J_5^^^+i/2). (32) 

We used a normalized grid parameterized by the 
fixed domain variable Z = z/h{t). This will make 
it possible to compare in a fixed frame the results at 
different times and depths and for different values 
of the dimensionless parameters. This transforma- 
tion maps the basement of the basin to Z = and 
the basin top to Z = 1. Numerical results are pre- 
sented and explained in the following section, and 
arc compared with the real data in the rest of the 
paper. 

5 Numerical Simulation and 
Comparison With Real 
Data 

The main aim of the numerical model is to show 
how the model equations behave, to test the va- 
lidity of the model, and to predict or simulate 
real world situations. In principle, a good model 
should be able to simulate many realistic features 
when its physical parameters are appropriately cho- 
sen, but in reality, the accuracy of the simulation 
is comprised by a lack of information about the 
real system. By changing the different parameters 
(A, A, (3, S, m, n and time t) as well as the boundary 
conditions {4>o and heat flux), we can in principle 
get a best fit to the real data although such a choice 
of such parameters is not unique because different 
combinations of parameters in the parameter space 
may correspond to fit the real data equally well 
with a same range of deviations. So the choice of 
the parameter is still rather arbitrary although the 
knowledge of the geology is used in attempt to get 
a more realistic grouping of papameters. 

The borehole log data (designed as LOG I) from 
the DQ-151 borehole in South China Sea Basin in 
China was used in the present simulations. DQ- 
151 has a depth of 6100 m and its porosity, volume 
fractions of smectite and illite, permeability and 
temperature distribution were surveyed at nearly 5 
meter intervals. The data used below, however, has 
been smoothed by averaging over intervals of 250 
m, so that only large scale features are revealed. 



However the profiles are still detailed enough to 

test the model presented in this paper. 

The results of the numerical simulations are com- 
pared with the data in Figure 1. The rescaled 
height Z = z/h(t) varies from to 1 and corre- 
sponds to a variation of 6100 m from the basement 
to the ocean floor. In simulating LOG I, the val- 
ues of the parameters which best fit the data were 
A = 250, S = 2.4, m = 8.3, n = 1.46, P = 2.4, 
t = 7.3, and ©c = 3.4. By using the typical length 
scale d = 970 m, ko « 1.22 x 10^^^ m^, we can es- 
timate the sedimentation rate at deposition to be 
rh, » 0.77 X 10"" m/s w 230 m Ma'^ The scal- 
ing time is about 3.3 Ma. and thus the dimensional 
time is about 24 Ma. For a fixed value of G = 1.0, 
the estimated viscosity is ^0 ~ 5x 10^^ N s m~^. We 
can see that the porosity near the top of the basin 
decreases nearly exponentially with depth. It can 
be approximately expressed as 

= 0oe-^»('"'-^\ (33) 

where Cg w 0.93 is the compaction index. Poros- 
ity reduces rapidly from 0.4 to 0.1. This expo- 
nential profile is consistent with the equilibrium 
profile given by Fowler and Yang (1998) found us- 
ing asymptotic analysis of the nonlinear reaction- 
diffusion equation for fast compaction in sedimen- 
tary basins. In the dynamic balance of sedimention 
and compaction, the porosity has reached equilib- 
rium in LOG I in the top 2100 m of the basin. 
Fowler and Yang (1998) also predicted that equi- 
librium even for a compacting basin may reach a 
depth of 



which is about 700m in the present case. The near 
equilibrium at the top means that (a) compaction 
in LOG I is mainly poroelastic and may reach the 
equilibrium state at depths much greater than that 
predicted by asymptotic analysis, and (b) no over- 
pressuring occurs in the shallow region above about 
2000 m, which is consistent with the drilling log 
data. However, in much deeper regions, the poros- 
ity profile is no longer exponential. A transition 
occurs at a depth of about 2400 m below which the 
porosity becomes nearly uniform. This means that 
compaction is in a non-equilibrium state and over- 
pressuring starts to build up. Compaction gradu- 
ally becomes viscous due to the mechanism of pres- 
sure solution which operates at higher pressures 
and temperatures than those prevailing at depths 
of less than 2000 m. 

Figure 2 shows the computed temperature pro- 
file and the corresponding data. The values used 
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in this figure are A = 150, 7 = 0.031 (degree/m) 
with the other values being the same as in Figure 
1. It is clearly seen that the temperature profile is 
essentially governed by conduction and evolves on 
a faster time scale than porosity reduction. This 
is because ~ —1 when A 1 and K 1. 

Consequently, the temperature distribution is vir- 
tually a linear profile versus depth. Thus we can 
take Q = {h — z). 

Figure 3 compares the computed volume frac- 
tions of different mineral species and with data 
from the borehole DQ-151. We can see that the 
volume fractions of smectite remains virtually con- 
stant from the ocean floor to a depth of about 2100 
m, and suddenly decrease at a depth of about 2400 
m, and becomes negligible at a depth of about 
5100 m. Assuming a thermal gradient of nearly 
30.7 C/km, the mineral reaction commenced at 
nearly 83°C, and it ended at nearly 160°C. As 
smectite starts to decrease, the volume fraction of 
illite starts to increase and it reaches its maximum 
at the aforementioned depth of about 5100 m. The 
reaction region has a thickness of about 2500 m and 
a temperature variation of nearly 77°C. The results 
in this paper are quite consistent with the general 
conclusions drawn from other data (Abercrombie 
et al, 1994) from oceanic and sedimentary basins. 

For a linear temperature distribution with depth 
(Yang, 2000b) has analysed the nonlinear equation 
asymptotically and concluded that the convection 
term d{(j)mU^) / dz <C 1, so that we have 



dt 



,l3{h-z-z')^ 



(35) 



and thus an approximate solution for the volume 
fraction of smectite is 

= <i)^o exp[-ie-'5(''-^-^*)], (36) 

where z* is the location of the centre of the reaction 
region, and h is the total depth of the basin. This 
solution is remarkably accurate in describing the 
features within the reaction window. The centre of 
the reaction window is 



(37) 



This clearly means that the higher the sedimenta- 
tion rate rhg, the higher the critical temperature 
T* of the mineral reaction, the deeper the reactive 
region z* , and vice versa. A change of 2 orders 
in sedimentation rate (m^) will cause a shift of 
by 2 (equivalently ~ 60° C) (with other parame- 
ters unchanged). In addition, the thickness of the 



reaction region is the order of (5//3)d. A typical 

value of /3 « 2.4 gives 2000 m (with d = 970 m), or 
equivalently a temperature range of ~ 63°C, which 
is quite close to the range observed in real data. 



6 Conclusion 

The present model of mechanical compaction, pres- 
sure solution, and the smectite to illite mineral re- 
action in hydrocarbon basins uses a rheological re- 
lation which incorporates both poroelastic and vis- 
cous effects and considers mineral reaction as a one- 
step dehydration mechanism in a 1-D compacting 
frame. The nondimensional model equations are 
mainly controlled by four parameters A, the viscous 
parameter S, the thermal conduction parameter A, 
and the reaction parameter TZ. The highly nonlin- 
ear coupled partial differential equations have been 
solved using a two stage implicit method which is 
quite robust for the present governing equations. 

The numerical simulations have shown that 
porosity-depth profile is near exponential at shal- 
low depths and then undergoes a transition to a 
nearly uniform porosity. This is because of the 
large exponent m in the permeability law k = 
(0/^0)™ J so that even if A 1, the product Xk may 
become small at sufficiently large depths so that 
the hydroconductivity becomes very small and thus 
compaction proceeds extremely slow (Fowler and 
Yang, 1998). Meanwhile, since the bulk viscosity 
depends on the porosity in a form of ^ = Co(?5'o/'/')" 
(typically n « 2), a decrease to sufficiently low 
porosity will give an increase to a very high vis- 
cosity. The transition from the higher porosity to 
a nearly minimum porosity is now associated with a 
sharp change of lower viscosity to a higher viscos- 
ity, and this implies a transition from poroelastic 
to viscous behavior. This is because of the Be- 
low this transition region, the porosity is usually 
uniformly small and viscous compaction essentially 
dominates. The sudden switch from poroelastic to 
viscous compaction is often associated with a sharp 
increase in pore pressure and a low permeability re- 
gion where the mineralized seal may be formed. As 
viscous compaction proceeds, porosity and perme- 
ability may become so small that fluid gets trapped 
below this region, and compaction virtually stops. 
Comparison with borehole data has shown very 
good agreement between computed porosity pro- 
files and smoothed borehole data. 

The temperature evolution and thermal history 
are mainly controlled by A, which is the ratio of 
the time scale of thermal conduction to the time 
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scale of the sedimentation rate. Observations sug- 
gest that the temperature profile is purely conduc- 
tive and evolves on a faster time scale than porosity 
reduction, that is to say, A 1 so that ^ » —7 
when K = const. This is consistent with the linear 
temperature profile generated by our model. 

The smectite to illite mineral reaction is charac- 
terised by the reaction parameter 7?., which may be 
defined in terms of a critical temperature 9c. This 
study reveals that mechanical compaction, which is 
controlled by the strata permeability and sedimen- 
tation rate, is the most important geological fac- 
tor in altering the porosity evolutions, while chem- 
ical compaction, controlled by mineral reaction and 
pressure solution, causes only small changes in the 
porosity. The first-order dcihydration model is a 
good approximation since it describes the extent 
of progress of the overall smectite-to-illite transfor- 
mation and still reproduces many essential features 
of the smectite-to-illite process if the appropriate 
reaction rate laws are used based on the known 
physics and chemistry from experimental studies. 
Naturally, more work is needed on more realistic 
formulation 2-D and 3-D compaction and mineral 
reactions in sedimentary environments. 
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Figure captions 

Figure 1. Comparison of computed porosity (solid 
curve) with borehole data (depicted as circles) 
taken from the South China Sea basin. Z = z/h{t) 
is the scaled height. Compaction in the top region 
of the basin is essentially at equilibrium (for 
A f« 250, S = 2.4, m = 8.3, n = 1.46, t = 7.3). 
The porosity profile decreases nearly exponentially. 

Figure 2 Temperature profile for A — 150 with all 
other values are the same as in Fig. 1. Numerical 
simulation (solid curve) is very consistent with the 
data from Log I. 

Figure 3 Comparison of numerical results (solid 
curves) with data from Log I. The volume fractions 
of smectite and illite change dramatically in a very 
thin reaction region. The values used in the simu- 
lations are A = 250, t = 7.3, 6^ = 3.4, ^ = 2.4. 
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